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X-ray Free Electron Lasers (XFEL) are revolutionary photons sources, whose 
ultrashort, brilliant pulses are expected to allow single molecule diffraction 
experiments providing structural information on the atomic length scale. This 
ultimate goal, however, is currently hampered by several challenging questions 
basically concerning sample damage, Coulomb explosion and the role of non- 
linearity. By employing an original ab-initio approach, as well as exceptional 
resources of parallel computing, we address these issues showing that accu- 
rate XFEL-based single molecule imaging will be only possible with ultrashort 
pulses of half of femtosecond, due to significant radiation damage and the for- 
mation of preferred multi-soliton clusters which reshape the overall electronic 
density of the molecular system at the femtosecond scale. 



Introduction Within the next few years three X-ray Free Electron Lasers (XFEL), namely the 
Linac Coherent Light Source (LCLS) at Stanford (1), the European XFEL in Hamburg (2) and 
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the Japanese XFEL at RIKEN (3), will start the operations. These facilities will generate the 
brightest X-ray pulses in the world, with atomic scale wavelength (> 0.05 nm) and maximum 
peak powers (0.1 — 0.5 TW) order of magnitudes beyond the current capabilities of synchrotron- 
based X-ray sources. Throughout history, the discovery of new light sources has given rise to 
decisive steps forward in science, and the new XFEL facilities hold the promise of making 
accessible a new fundamental physical domain, thus establishing a new era in science and re- 
search. This has led to the flourishing of an incredibly large interest around the XFEL in recent 
years (4-15). Among the various challenges proposed — which range from extreme investiga- 
tions in high-energy physics to ultrafast chemical analysis — the possibility to perform coherent 
diffraction imaging at the atomic scale has stirred particular interest (4, 7,9, 10, 12, 13, 16, 17). 
The difficulties in determining the structure of proteins that cannot be crystallized, in fact, is 
one of the leading challenges for structural biology today, and the extremely brilliant, ultra- 
short, bursts of the XFEL have the potential to get around this problem; on the one hand the 
XFEL pulse duration (femtoseconds), is the natural time scale of elementary atomic interac- 
tions and it is expected to open the study of the evolution of fundamental processes, which 
could otherwise be studied only indirectly; on the other the XFEL could have the power to 
enable diffraction imaging of single molecules, thus overcoming conventional crystallography 
that works essentially through linear amplification (by observing discrete Bragg peaks resulting 
from the interference of regular arrangement of atoms). The initial steps for this new frontier of 
research have begun at the FLASH facility (18), where imaging experiments have successfully 
captured a time-series of snapshots showing the nonequilibrium, ultrafast transient dynamics of 
a nanostructured solid: the laser ablation of a silicon film (19) — with a spatial resolution of 
50 nm and a temporal resolution of 10 ps — and the explosion dynamics of polystyrene spheres 
with diameter of 140 nm (10). However, to reach the ultimate goal of single molecule imaging, 
things turn out to be more complicated. Aside from the most fundamental questions concerning 
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the phenomenology arising from high intensity photon beam-matter interactions (which are not 
fully exploited) there is a practical aspect that needs to be taken into account. Atomic scale 
imaging of particles, in fact, will be lensless, with the far field diffraction pattern scattered by 
the molecule first recorded on a camera and then analyzed through a computer algorithm, which 
reconstructs the molecular image (4, 13, 16, 17, 20-22). To theoretically support this approach, 
several empirical models of ultrashort photon pulse interactions with matter have been pursued 
previously (23-27). However, owing to the large number of approximations employed, none 
have given a definitive answer to the basic question: is it possible to collect the scattering signal 
from an electronic ground state, to permit the reconstruction of the molecular structure? Within 
all these models, in fact, the primary interest is in the dynamics of atomic nuclei, while elec- 
trons are treated with very crude assumptions, ranging from simple rate equations to scattering 
cross sections, all of them derived in the framework of perturbative analysis of noninteracting 
electrons subjected to low external excitations. XFEL imaging, on the contrary, relies on the in- 
terplay between intense X-rays pulses and electron wavefunctions that, in turn, interact through 
nonlinear and extremely nonlocal (in space) functionals of the overall electron density (28) . An 
additional challenge to existing approaches also comes from the latest results obtained at LCLS, 
where ultrashort pulses less than 10 fs length have been generated. Simulation indicates these 
pulses may be as short as 2 fs with peak power levels of perhaps 500 GW. (29). The simulation 
provides novel insights in the internal structure of XFEL radiation, unveiling the presence of 
trains of incoherent energy bursts with mean time duration as short as 200 as, whose interaction 
with matter cannot be modeled without a first-principle theory. 

The frontier of atomic resolution imaging is still largely debated with the most problematic 
issues unsolved. Open questions include the photoionization and subsequent Coulomb explo- 
sion dynamics of a single molecule subjected to intense XFEL radiation (and in particular, if 
it could survive long enough to allow the observation of a useful diffraction signal) (9), what 
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is the sample far field scattering pattern and if it encodes interference fringes with sufficiently 
high contrast to infer structural information on the material and, last but not least, what is the 
role of nonlinearity that, on the basis of recent theoretical calculations made on plasma models, 
is expected to affect the dynamics of high energy XFEL pulses (30). 

The aim of this work is to answer to these questions. On one hand, to define the basis for present 
and future research in XFEL science, we develop an ab-initio model describing the nonlinear 
interaction of intense XFEL beams with matter; on the other, to address the above-mentioned 
open questions we employ a numerical parallel approach with the highly scalable code GZilla 
(especially written and optimized for this problem) performing simulations with a variety of 
atoms and simple molecules. Our results are twofold: 

i) we addressed the problem of sample radiation damage and demonstrate that the full pho- 
toionization of simple molecules occurs in a few (< 5) of femtoseconds, with valence electrons 
escaping in 200 — 300 attoseconds; by a further comparison with the photoionization time of 
the hydrogen atom, we elucidate the fundamental role of nonlinear electron interactions in the 
photoionization dynamics. 

ii) we investigated the coherent imaging capabilities of XFEL sources, collecting snapshots of 
integrated far field diffraction patterns (as it would be retrieved by a standard camera), ions 
positions, electrons and electromagnetic energy density of molecules illuminated by ultrashort 
pulses (both in the femtosecond and attosecond regimes). We highlighted the existence of com- 
petitive dynamics on the femtosecond time scale, sustained by the strongly nonlinear nature of 
the electrons interactions, whose effect is to appreciably alter the electron density of the molec- 
ular system. 

The main conclusion of this work is that, albeit some structural information information is still 
retrievable in the femtosecond domain, accurate single molecule imaging will be only possible 
in the sub-femtosecond regime. 
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A first principle model and a state-of-the-art parallel solver To derive a first principle 
model of an ensemble of nonrelativistic atoms under the presence of a time-dependent electro- 
magnetic field, we resort to the following quantum Hamiltonian H (31): 
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with atomic nuclei (electrons) defined by the charge density operator n n (n e ), the momentum 
Pn (PeX the atomic number Z (elementary charge e) and mass m n (m), while the electromag- 
netic (e.m.) field by the electric E and magnetic H field (or equivalently by the vector potential 
A), with e and fi being dielectric and magnetic constants, respectively. Equations of motion 
for both matter and e.m. field can be found directly from (OQ) through the application of the 
Poisson bracket operator (supporting online material). At variance with previous theoretical 
and numerical studies, here all the quantities appearing in the quantum Hamiltonian (i.e., nu- 
clei, electrons and e.m. field) are treated as dynamical variables, evolving through the nonlinear 
set of equations arising from the Hamiltonian (OQ) itself (supporting online text). 
The numerical simulation of the dynamics arising from (OQ) is accomplished with the paral- 
lel code GZilla (supporting text online), which combines classical Density Functional Theory 
(DFT) (32) for the computation of the system ground state, and TDDFT coupled to Finite- 
Difference Time-Domain (FDTD) (33) and Molecular Dynamics (MD) codes for the time- 
dependent analysis of (OQ). GZilla has been parallelized with the Message-Passing-Interface 
(MPI) standard, which allows top performance in processor scalability, tested and verified on 
up to 4096 processors to date. The results presented have been derived by 15 million of single 
CPU computational hours, representing one of the longest simulations ever reported; thanks 
to the code scalability and to the availability of parallel computing resources, more than 1700 
years of single cpu computations have been completed in less than 3 months. 
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Setting the ground state In order to provide an adequate and experimentally interesting 
selection of molecules, we look at the four commonly-found elements of organic chemistry 
(H, O, N, C) and study their stable configurations, ranging from the most elementary com- 
pound to the simplest combination involving all of them (Fig. [T]). Ordered by an increasing 
degree of electronic complexity, we considered the hydrogen atom H (Fig. CD-A), water H 2 
(Fig. ffl-B), methane CH A (Fig. \T\-C), nitrogen (Fig. ED-D) and finally the isocyanic acid 
HNCO (Fig. CD-E) (34). The chosen molecules exhibit a mean bond length of about 0.1 nm, 
which matches well with the XFEL shortest wavelength 0.1 nm), and different sizes of exter- 
nal electron clouds, allowing the investigation of the interplay between nonlinearity and sample 
geometry at various spatial scales. 

Molecules in high intensity e.m. field: radiation damage and photoemission rates The 

first series of ab-initio simulations are devoted to the study radiation damage of samples sub- 
jected to intense XFEL radiation, which is a key issue to calibrate XFEL experiments. The spa- 
tial distribution of the input source has been modeled with a Gaussian spatial profile, bounded 
to 100 nm, a focusing which appears feasible with the existing technology. To investigate the 
interplay between nonlinearity and geometry, we first considered a continuous (cw) excitation. 
Parameters for wavelength and power have been chosen to match the XFEL in the high energy 
regime; in particular, for wavelength A well below 1 nm, we considered the SASE 3 configu- 
ration of the European XFEL (A = 0.4 nm, power P = 150 GW) (2). Figure (0Q) shows the 
photoionization of H, H 2 0, CH4, N 2 and HNCO, displaying the time evolution of the to- 
tal number of electrons in the computational box (a cube of 1.5 nm side for iV 2 and 1 nm for 
all the others). As seen in Fig. [Q the sample radiation damage is extremely high: with the 
only exception the H atom, large ionization occurs in a few of femtoseconds with the HNCO 
molecule exhibiting a radiation damage (percentage of electrons lost) of 98% at a time t = 4 
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fs, showing two different characteristic exponential ionization rates versus time. As a result, at 
the single femtosecond scale (much before a measurable Coulomb explosion of the nuclei), the 
dynamics is completely dominated by electrons shaking up; since the latter are the only source 
of scattered photons (nuclei, in fact, are too heavy to produce any appreciable scattering), any 
imaging-oriented XFEL application should be realized within ultrashort (< 1 fs) pulses. 
A second key observation comes from the different behavior of H, which indicates that non- 
linearity, rather than geometry, outweigh the photoionization process. In fact, as displayed in 
Fig. |2K, the internal electrons of H 2 and HNCO, whose wavefunctions have spatial extent 
comparable to those of H, escaped much before those of H showing that electron nonlinear 
interactions actually dominate the system evolution. A physical interpretation comes from the 
theory of solitons (35); the electron ground- state, in fact, being a solution of the dynamical 
equations arising from (OQ) with constant amplitude and linear phase evolution, can be regarded 
as a a set of M solitary waves interacting through a functional of the electron density (i.e., the 
V xc potential), which plays the role of a nonlinear Kerr-like response. Such nonlinearity is well 
known to foster processes with a considerable amount of radiation emission, as soon as solitary 
waves are perturbed from their ground state and forced to overlap. Therefore, when the XFEL 
pulse actually breaks in and upsets the atomic system, it initiates a process of radiation emission 
that gets amplified by the system nonlinear response through the interaction among electronic 
solitary waves. It is clear that this mechanism does not play a role for hydrogen, which has just 
one electron, resulting in a much lower photoionization rate with respect to heavier atoms. 

Single molecule imaging: a benchmark with the isocyanic acid HNCO The most challeng- 
ing (and unsolved) issue concerns the possibility of recording the far field molecular interference 
image on a camera that, can not employ any time gating and will detect scattered photons within 
the duration of the X-ray pulse. To answer to this question, we performed a series of simula- 
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tions collecting the time domain far-field pattern (33), integrated in time as it would be done by a 
standard CCD, scattered by an HNCO molecule irradiated with an high energy XFEL pulse of 
Gaussian spatial distribution (100 nm spot-size) under different ultrashort regimes. We begin by 
considering the XFEL short pulse operation, recently characterized at the LCLS facility, which 
is expected to produce 2 fs pulses with peak powers of P = 0.5 TW at A = 1.5 nm (29). This 
regime allows a direct comparison with the XFEL technology accessible in the near future. The 
set of Figs. dU)-© provides a comprehensive analysis by showing the photoionization (Fig. [3] 
A), the distance between the nuclei (Fig. [3]B), the electron/nuclei dynamics and the angular far 
field (Movie SI and movie snapshots Fig. H]) time evolution of an HNCO molecule irradiated 
by a XFEL short pulse, constituted by a random series of Gaussian- shaped incoherent energy 
bursts (single burst length m 200 as, overall length 2 fs), with peak power P = 500 GW and 
wavelength A = 0.15 nm (Fig. |3]A dashed line). As seen in Fig. |3]A, each XFEL energy burst 
pushes the HNCO ionization a step forward, appreciably damaging the molecule much before 
the arrival of the main pulse peak (Fig. [3] A, dashed line). At the point of maximum XFEL 
intensity 1 fs), in fact, the HNCO molecule has lost three electrons; after the pulse peak, 
HNCO ionization proceeds at a slower rate and at the end the radiation damage reaches a value 
of 23%, with more than five electrons escaped from the molecule. Figure [3] B illustrates the 
consequences of molecular damage on nuclei dynamics. In particular, the hydrogen atom is no 
longer bound to the structure (after 400 as) and initiates the process of Coulomb explosion. The 
nuclei pairs C — N and C — O, conversely, provide opposite dynamics: while the latter exhibit 
a tendency towards explosion, the former are not affected by this process and evolve with a con- 
stant bond length within the XFEL time window. In the framework of the previous nonlinear 
wave interpretation of the dynamics, such a different evolution is the hallmark of the generation 
of a robust multi-soliton cluster, with a consequent spatial reshaping of the electron density due 
to the nonlinear mixing of C and iV electron waves. Movie SI and Figs. |4]A-D provide a 
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clear illustration of this nonlinear process. At the beginning (t=0 fs), the carbon atom shows 
no preferential coupling to either nitrogen or oxygen, as witnessed by the symmetric density 
distribution which equally affects C and N neighbors. However, as soon as the hydrogen atom 
looses its coupling to the structure (Fig. |4]C), the situation dramatically changes: the carbon 
and nitrogen atoms interact and evolve towards a robust multi-soliton cluster, as demonstrated 
by the strongly asymmetric density distribution which survives at the end of the XFEL pulse 
(Fig. H]D). 

In conclusion, even by the use of short pulses and in the presence of negligible Coulomb ex- 
plosion, the electronic radiation damage of molecular samples is still an issue: on one hand, it 
deprives the molecules of an appreciable number of electrons, on the other, it initiates a change 
of the electron density through the formation of preferred multi-soliton clusters, which tend to 
oppose sample ionization. The net effect of these two different dynamics (both sustained by 
the strongly nonlinear nature of the light-matter interaction process) is a qualitative and quan- 
titative reshaping of the electron density which, in turn, leads to the inability to reconstruct the 
molecular structure from the far field diffraction pattern. In other words, although the few fem- 
tosecond time scale is short enough to avoid a significant Coulombian explosion of the nuclei, 
the sub-femtosecond (« 0.5 fs) electron density reshaping is still an issue for single molecule 
imaging experiments. 

A general improvement is possible by moving to the attosecond domain, as expected on the 
basis of Fig. [3] A, which predicts considerably smaller radiation damage after the first energy 
burst (f=s 0.5 fs). Figure [5] summarizes the results of attosecond molecular diffraction: by em- 
ploying a 400 as long ultrashort XFEL pulse, the far field diffraction pattern (despite the short 
illumination) shows good fringe contrast (Fig. [5] A), much superior to the pattern obtained in 
the fs regime (Fig. |4]D), minimal radiation damage (8.5% as calculated from Fig. 03) and no 
tendency toward nuclei explosion (Fig. [5]C). Figure [5)3, showing the difference between the 
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previous far field diffraction pattern and the HNCO diffraction from the ground state (i.e., in 
the absence of any molecular damage), highlights the minimal effects of radiation damage in 
the attosecond diffraction regime (Fig. [5]A,D). The sub-femtosecond domain, in other words, is 
sufficiently fast to minimize the nonlinear interaction dynamics of electron matter waves (Fig. 
[5K) and to maintain all the nuclei bonded together (Fig. [5]C). A computer reconstructed image 
from an attosecond diffraction pattern would contain all the information of the original electron 
density and will therefore be a true representation of the original sample. 
In summary, in both short and ultrashort regimes the Coulomb explosion of the nuclei is not an 
issue for single molecule imaging, but the radiation damage (due to the high energy carried by 
the XFEL pulse) do offer some limitations to Angstrom-scale imaging. In particular, competi- 
tive dynamics originated and sustained by nonlinear electron interactions lead to a distortion of 
the electron density distribution in the XFEL femtosecond regime. Albeit some structural infor- 
mation is still retrievable by employing femtosecond XFEL pulses, a computer reconstructed 
image will no longer be a correct reproduction of the original sample due to the unbinding of 
some molecular nuclei and the formation of preferred multi-soliton clusters. Such dynamics 
are avoided in the attosecond domain, which is sufficiently faster then the characteristic scale 
of such effects, preserving the original electron distribution and the resulting diffraction pattern 
from loosing any important information. 

Dream or reality? With the opening of LCLS and the collection of its initial experimental 
results, the XFEL has become a reality of unprecedented importance. In this paper, we have 
derived an original and accurate description of XFEL beam interaction with matter, able to 
provide ab-initio results and expected to settle the basis for present and future theoretical in- 
vestigations in XFEL science. When applied to the one of the most promising applications of 
XFEL sources, our results predict that some work still needs to be done in order to reach the goal 
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of single molecule imaging and Angstrom scale microscopy, but we foresee that such frontiers 
may be just behind corner, especially in conjunction with the fast theoretical and experimental 
advances in the field of attosecond physics (15). 

References and Notes 



1 . http://lcls.slac.stanford.edu/ 



2. http://www.xfel.eu/ 



3. http://www.riken.jp/xfel/eng/index.html. 

4. K. J. Gaffney, H. N. Chapman, Science 316, 1444 (2007). 

5. T. Shintake, et al, Nature Photon. 2, 555 (2008). 

6. M. Dunne, Nature Physics 2, 2 (2006). 

7. J. D. Brock, Science 315, 609 (2007). 

8. D. Clery, Science 322, 1769 (2008). 

9. H. N. Chapman, Nature Materials 8, 299 (2009). 

10. H. N. Chapman, et al, Nature 448, 676 (2007). 

11. S. P. Hau-Riege, et al, Phys. Rev. Lett. 98, 145502 (2007). 

12. A. Cavalleri, Nature 448, 651 (2007). 

13. J. Kirz, Nature Physics 2, 799 (2006). 

14. P. Emma, et al, Phys. Rev. Lett. 92, 074801 (2004). 

11 



15. F. Krausz, M. Ivmo\,Rev. Mod. Phys. 81, 163 (2009). 

16. K. A. Nugent, Nature Physics 5, 17 (2009). 

17. R. Fung, V. Shneerson, D. K. Saldin, A. Ourmazd, Nature Physics 5, 64 (2009). 



18. http://flash.desy.de 



19. A. Barty, Nature Photonics 2, 415 (2008). 

20. C. G. Schroer, et al, Phys. Rev. Lett. 101, 090801 (2008). 

21. J. Miao, P. Charalambous, J. Kirz, D. Sayre, Nature 400, 342 (1999). 

22. S. Eisebitt, et al, Nature 432, 885 (2004). 

23. Z. Jurek, G. Faigel, M. Tegz, Eur. Phys. J. D 29, 217 (2004). 

24. R. Neutze, R. Wouts, D. van der Spoel, E. Weckert, J. Hajdu, Nature 406, 752 (2000). 

25. S. P. Hau-Riege, R. A. London, A. Szoke, Phys. Rev. E 69, 051906 (2004). 

26. S. P. Hau-Riege, H. N. Chapman, Phys. Rev. E 11, 041902 (2008). 

27. M. Bergh, G. Huldt, N. Timneanu, R. N. C. Maia, Filipe, J. Hajdu, Q. Rev. Biophys. 41, 
181 (2008). 

28. E. Runge, E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). 

29. J. Frisch, Short pulse operation, Science Drivers for Hard X-Ray Upgrades to LCLS Work- 
shop, Menlo Park, California, July 29-31 2009. 

30. A. Fratalocchi, C. Conti, G. Ruocco, F. Sette, Phys. Rev. B 11, 245132 (2008). 

31. R. Loudon, The Quantum Theory of Light (Clarendon Press, Oxford, 1983). 

12 



32. R. M. Martin, Electronic Structures (Cambridge Press, Cambridge, 2004). 

33. A. Tafiove, S. C. Hagness, Computational Electrodynamics: the finite-difference time- 
domain method (Artech House, 2000), third edn. 

34. Electron densities are rendered with the volume mapping representation of the visualization 
toolkit (vtk) library (http://www.vtk.org), which maps each intensity point within the three 
dimensional cartesian grid to a specific color with a particular transparency, projecting the 
resulting image to a viewpoint provided by a camera observer. 

35. Y. S. Kivshar, G. P. Agrawal, Optical solitons (Academic Press, New York, 2003). 

36. M. A. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, Time-Dependent Density Functional 
Theory (Springer, Berlin, 2006). 

37. F. Schwabl, R. Hilton, A. Lahee, Advanced Quantum Mechanics (Springer, Berlin, 2008). 

38. This research is done in the framework of The SolarPaint project (www.solarpaintproject.org), 
an activity funded by Award No. KUK-F1-024-21 (2009/2012) made by King Abdullah Uni- 
versity of Science and Technology (KAUST). We are indebted to KAUST (Aron Ahmadia, 
Moody Altamimi, Bruce Guile, Dinesh Kaushik, David Keyes, Laura Meany, Richard Orme) 
and KAUST-IBM partnership (Chris Beyer, Pavan Kumar, Cindy Mestad, Fred Mintzer, 
Dave Singer, Julie Sherlock) for the large parallel resources provided through KAUST Sha- 
heen, IBM WatsonShaheen and IBM Rochester Blue-Gene supercomputers. Code devel- 
opment has been done at the Max-Plank IBM Power6 system, thanks to the Distributed 
European Infrastructure for Supercomputing Applications (DEIS A) European award given 
to the Ab-initio Coulomb Explosion Simulations at X-rays (ACES-X) project. We are very 
grateful to the CINECA team (Cristiano Calonaci, Andrew Emerson, Giovanni Erbacci) for 
the parallel support in the enabling stage of ACES-X. We finally acknowledge illuminating 

13 



discussions with Massimo Altarelli (DESY European XFEL) and Jerry Hastings (Stanford 
LCLS) concerning the practical capabilities of XFEL light sources. 

Supporting Material 

SOM text 
Movie SI 

References S2, S3, S4 



14 



Figure 1: Three dimensional images of ground state molecules, with nuclei displayed as rigid 
spheres of arbitrary size, of (A) H, (B) H 2 0, (C) CH 4 , (D) N 2 and (E) HNCO (the bottom-left 
part of each panel shows the corresponding number of electrons). 
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Figure 2: (A) Number of 'bounded' electrons versus time for H, H 2 0, CH 4 , N 2 , and HNCO 
molecules illuminated by a 100 nm-waist Gaussian beam of cw XFEL radiation with power 
P = 150 GW and wavelength A = 0.4 nm. 
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Figure 3: Time evolution of electrons number (A, solid lines) and nuclei distance (B) of 
HNCO irradiated by a random train of incoherent energy bursts of peak power P = 500 
GW, wavelength A = 0.15 nm and overall time length T=2 fs (A, dashed line). 




Figure 4: (A)-(D) Far-field scattered angular pattern (red to yellow colormap), nuclei position 
and electron density (blu to yellow colormap) time evolution of an HNCO molecule irradiated 
by a short XFEL pulse of 2 fs length, wavelength A = 0.15 nm and peak power P = 500 GW 
(Fig. [3] A dashed line). 
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Figure 5: Attosecond diffraction results: (A) Angular Far-field, nuclei position/electron density, 
(B) photoionization and (C) nuclei distance evolution of an HNCO molecule irradiated by a 
400 as XFEL pulse of wavelength A = 0.15 nm and peak power P = 500 GW (B, dashed line); 
(D) Image difference between the far field of (A) and the one generated by the ground state of 
HNCO in the absence of any radiation damage. 
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Supporting Material 

1 Supporting Text 
1.1 Equations of motion 

Atomic nuclei, due to their extremely small size (~ 10~ 15 m) with respect to hard X-ray wave- 
lengths (> lCT 10 m), are seen as point objects from the XFEL, hence, the quantum expecta- 
tion of the nuclei density (\l/|n n (r)|\l/) = J2i$( r — R-i) (averaged over the many body state 
\^) = IV^Ix) composed by electrons \ip) and nuclei \x)) factorizes into a series of Dirac delta 
centered at the nuclei positions Rj. The latter, and their conjugate momenta Pj, evolve accord- 
ing to the following Poisson brackets: 

r)Ti 

^ = {R*<*|H|*>}, 

-^l = {P i ,<*|H|«r>}, (2) 

with i = 1, N. After some algebra, we end with [J: 

<9 2 R 

Mi-^ = Z ie [(E Ti + v, x Bi) - V Ri 0(R)], (3) 
with E Ti = E T (Rj) [the subscript T denotes the transverse part] and the potential 4>: 

n l) 4 7 re V^|R i -R i | J \Ri-r\J' 

being Zi the atomic number of the i— th nucleus and n the quantum expectation of the electron 
density. 

The evolution of the electron many body state in the Schrodinger representation, is de- 
scribed by: 

in^lV) = H|v>) (5) 



'Throughout this article, we assume to work in the Coulomb gauge. 
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The latter, in its direct form, is unsuitable for a numerical integration due to a dimensionality 
bottleneck problem. The memory workload for storing the wavefunction ip, in fact, grows 
exponentially as L M for an M— body system with size L. By using a tiny grid of 128 points in 
each spatial dimension (typical simulations fall in the range 256 — 1024 points), more than 60 
millions of terabytes would be required for storing the wavefunction of just 3 electrons, which 
makes the problem unmanageable even in this simple situation. A more convenient formulation 
is provided by Time Dependent Density Functional Theory (TDDFT) (32, 36), which rewrites 
the single many-body problem © into a series of one-body problems coupled via a nonlinear 
potential V xc (the exchange-correlation potential), function of the electron density n, which 
describes the overall interaction among electrons. This formulation scales as L ■ M, and can 
be handled in a sufficiently-large supercomputing environment. As in Hartree-Fock theory and 
noninteracting quantum many-body dynamics (32, 37), TDDFT writes down the many body 
state as a Slater determinant: 

Wri,...,r M )) = ^ (6) 

with Sij = \ipj{vi)), and ipj = (v\tpj) playing the role of a single orbital whose evolution, 
according to Eq. dTJ, is found to be: 



Ti 2 „o e 2 / r , , n(r') ^ Z 4 \ _ e 2 . _ e 



— V 2 +tM dr'j^--^ - 1 , )+V xc +— A 2 -zft-A-V 
2m 47re n \J r — r , R,- — r / 2m m 



*l>j, (V) 



dt 

being A the electromagnetic vector potential and j running from 1 to M. 

The equations of motion of the electromagnetic fields, finally, are obtained from: 

^ = {q, (*|H|*», 

^ = {p, (*|H|«r», (8) 
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with p = — e E and q = A acting as momentum and position, respectively. Equations © yield 
the following nonlinear Maxwell's set: 



V x H + J 



dE _ 1 

dt e 
OH _ V x E 

dt [M 

V ■ E r = V ■ B = V • A = 0. (9) 



with a source J, representing the classical current of electrons and nuclei: 



N ZP 1 eA ' 
V -^^<5(r - Rj) (^|n e p e |^) n 

i=i 1 



(10) 



The set of equations represented by ©-©, © and (l9l)-(fT0l) forms the complete theoretical 
model of XFEL interaction with matter. The first principle XFEL model integrates molecular 
dynamics [Eqs. ©-©], Maxwell's [Eqs. (l9l)-(fT0l)l and hydrodinamics-like [Eqs. ©] equations 
in nonlinearly coupled platforms, which evolve on comparable time scales. A further compli- 
cation in the simulation of the set ©-([4]), © and (l9l)-(fT0l) originates from the mixing of both 
transverse and global quantities in Maxwell's equations (l9~l)-(fT0l). as well as the extremely small 
contribution arising from the scattered current J, which generate radiation fields several orders 
of magnitude lower than the XFEL beam. More specifically, according to the Helmholtz theo- 
rem (S2), the calculation of the transverse part of a generic vector field requires the evaluation 
of a 3D vector integral (which is very time consuming), while the extrapolation of such a small 
scattered field poses challenges even to the advanced total field/scattered field (TFSF) formula- 
tion for solving Maxwell's equations (able to isolate the scattered field contribution stemming 
from a single numerical grid). We addressed these problems by defining two different sets 
of Maxwell's equations, one dealing with the transverse field and one with the scattered con- 
tributions, exploiting the division naturally arising from the use of the Coulomb gauge (S2). 
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In particular, we considered the following equations modeling the evolution of the transverse 
beam: 

<9E T 1 

— — = -V x H T , 

dt e 

dH T _ V x E T 

dt fi 
dA 

~m ~ ~ Et ' 

V • E T = V • B T = V • A = 0. (11) 
and the next set for the dynamics of the scattered field: 



at e L 

<9H V x E 



(12) 



dt fi 

which is triggered by the current J generated by electrons and nuclei through (fTOl ) and depend- 
ing on the vector potential A of (PTTT) . Such a transverse field/scattered field formulation is 
particularly convenient since it allows to safely neglect the transverse scattered contributions in 
(ITT]) , as they are orders of magnitude lower than the ultraintense XFEL beam, yielding a set of 
equations that can be integrated very efficiently (no mixing with transverse and total fields) and 
which provide an accurate description of the scattered field arising from the molecules through 

CE2D. 

1.2 GZilla code implementation details 

GZilla features one of the most comprehensive simulator ever realized, composed to date of 
more than 120000 lines of parallel code. Its analysis capabilities range from DFT ground state, 
to time-dependent quantum molecular dynamics in the presence of electromagnetic fields of 
arbitrary form. Time-dependent far-field pattern analysis is also provided. GZilla implemen- 
tation has been developed by taking into account state-of-the-art methods and algorithms. In 
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particular: 

Maxwell's equations (l9l)-(fT0l) are solved by the FDTD technique within a Cartesian Yee grid 
with time marching Leapfrogging, Uniaxial Perfectly Matched Layers (UPML) and Total Field 
Scattered Field (TFSF) source formulation. The far-field pattern is calculated by the time- 
domain technique (33). Code parallelization is provided through spatial and time domain de- 
compositions. 

The time evolution of Schrodinger equations © is conversely committed to an original, un- 
conditionally stable Crank Nicholson propagator designed with graph theoretical approaches 
to minimize communications among processors and an original theory of Perfectly Matched 
Layers (PML) to absorb outgoing electron waves, combined with a Galerkin-based geometric 
Multigrid (S3) for the evaluation of the density integral appearing on the r.h.s of Eq. ©. Code 
parallelization is provided through a spatial domain decomposition strategy. 
Equations ©-©, finally, are solved by a velocity-verlet time stepping algorithm (32), paral- 
lelized by an original communicator splitting technique. 

The ground state DFT analysis is derived through a Rayleigh Quotient Multigrid (RQMG) 
eigenvalue solver (S4), within a Self-Consistent Field (SCF) iteration procedure (32) and an 
original guaranteed-reduction charge-mixing scheme. 
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